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Abstract. 


We have analyzed three years of radio tracking data from the MESSEN- 
GER spacecraft in orbit around Mercury and determined the gravity held, 
planetary orientation, and ephemeris of the innermost planet. With improve- 
ments in spatial coverage, force modeling, and data weighting, we rehned an 
earlier global gravity held both in quality and resolution, and we present here 
a spherical harmonic solution to degree and order 50. In this held, termed 
HgM005, uncertainties in low-degree coefficients are reduced by an order of 
magnitude relative to the earlier global held, and we obtained a preliminary 
value of the tidal Love number k 2 of 0.451±0.014. We also estimated Mer- 
cury’s pole position, and we obtained an obliquity value of 2.06 ± 0.16 ar- 
cmin, in good agreement with analysis of Earth-based radar observations. 
From our updated rotation period (58.646146 ± 0.000011 days) and Mer- 
cury ephemeris, we verihed experimentally the planet’s 3; 2 spin-orbit res- 
onance to greater accuracy than previously possible. We present a detailed 
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analysis of the HgMOOS covariance matrix, and we describe some near-circular 
frozen orbits around Mercury that could be advantageous for future explo- 
ration. 
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1. Introduction 


The gravity field of a planet provides fundamental information on the structure and 
evolution of the planet’s interior. In this paper, we report a new solution for the gravity 
field of Mercury, here termed HgM005, developed at the NASA Goddard Space Flight 
Center (GSFC) from nearly three years of radiometric tracking data from the MErcury 
Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER) spacecraft. 
The analysis presented here includes orbital observations between March 2011 and Eebru- 
ary 2014 spanning more than six Mercury years, or about 20 Mercury spin periods, as well 
as data from the three Mercury flybys in 2008 and 2009. In comparison to the previous 
MESSENGER models of Mercury’s gravity field [Smith et al, 2012; Genova et al, 2013], 
HgM005 is of higher resolution, to spherical harmonic degree and order 50, and includes 
refined force modeling that improves confidence in the solution. 

1.1. History 

Mercury, the innermost planet of the solar system, has been observed by humans since 
earliest time. Its proximity to the Sun made it challenging to study, and the planet 
remained largely a mystery even as our knowledge of other celestial bodies increased 
through ground-based observations and early spacecraft exploration. Mercury’s orbit 
around the Sun was characterized early, and the precession of its perihelion provided 
Einstein with an early experimental confirmation of general relativity, but the rotation 
period of Mercury was not reliably measured until 1965. Indeed, Mercury was earlier 
thought to be in a synchronous spin-orbit resonance, with a rotation period equal to 
its orbital period of ~ 88 days. Pettengill and Dyce [1965] reported a much shorter 
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rotation period, 59 ± 5 days, from ground-based radar measurements. Additional radar 
observations, and surface images by NASA’s Mariner 10 spacecraft [Klaasen, 1976] further 
confirmed that Mercury is in a 3: 2 resonance. 

The two equatorial flybys and one high-inclination flyby by Mariner 10 in 1974 and 
1975 also provided the hrst measurements of Mercury’s gravity held. Anderson et al. 
[1987] estimated the gravitational parameter [GM = (2.203209 ± 0.000091) x 10^^ m^ 
s“^, where M is Mercury’s mass and G is the gravitational constant] and quadrupole held 
[represented by the two second-degree terms in the spherical harmonic expansion of the 
gravitational potential, C 20 = (—2.68 ± 0.9) x 10“® and C 22 = (1.58 ± 0.8) x 10“®, dvr- 
normalized]. Unfortunately, the historical radio tracking data from the Mariner 10 hybys 
are no longer available because of a non-standard format and missing documentation. 

The MESSENGER spacecraft made three equatorial hybys of Mercury at an altitude 
of closest approach of ~ 200 km in 2008 and 2009. As reported by Smith et al. [2010], 
this geometry enabled good recovery of the difference in the equatorial moments of inertia 
[C 22 = (1-26 ± 0.12) X 10“®]. In contrast, the geometry was not favorable to the recovery 
of the gravitational polar battening {C 20 ), and Smith et al. [2010] cautioned that the 
estimated value of (—0.86 ± 0.30) x 10“® for that quantity was implausible given its 
implications for interior models. 

With ground-based radar observations obtained during 2002-2006, Margot et al. [2007] 
measured precisely the spin-axis orientation and the amplitude of the forced libration. 
These geophysical parameters provide important constraints on the interior structure of 
Mercury. The measured obliquity (2.11 ± 0.1 arcmin) from the spin-axis orientation 
provided observational evidence that Mercury is in or very near a Cassini state [Peale 
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et ai, 2002]. The large measured amplitude of the libration (~ 425 m at the equator) 
[Margot et ai, 2007; Margot, 2009] provided evidence that the outer solid shell of Mercury 
is decoupled from a liquid outer core. 

After MESSENGER was inserted into orbit about Mercury in March 2011, radio track- 
ing data led to a much improved determination of the low-degree held and the recovery 
of spatially resolved mass anomalies. From an analysis of the hrst several months of or- 
bital data, Smith et al. [2012] obtained a gravity held solution to harmonic degree and 
order 20, which they termed HgM002. Despite strong correlations between zonal har- 
monics because of MESSENGER’S eccentric orbit, the C 20 value was tightly constrained, 
to (—2.25 ± 0.01) X 10“®. C 22 was nearly unchanged from the hyby determination, at 
(1.25 ± 0.01) X 10“®. Modeling of interior structure was constrained by these updated val- 
ues and their improved uncertainties, and led to notable changes in the understanding of 
Mercury’s interior [Smith et ai, 2012]. The core fraction was revised still farther upward, 
and a solid FeS layer at the top of the liquid core was suggested as a means to reconcile 
the thin silicate shell with the large moment of inertia of the outer solid shell (to which 
the FeS layer would contribute). A more detailed analysis of the interior structure [Hauck 
et ai, 2013] with an updated value for the obliquity of 2.04 ± 0.08 arcmin [Margot et al, 
2012] showed that the FeS layer, although still allowed, is not required (see also Rivoldini 
and Van Hoolst [2013]). 

1.2. The MESSENGER Mission 

The MESSENGER spacecraft was launched on 3 August 2004 from Cape Canaveral. 
After an interplanetary trajectory more than six years long that included one flyby of 
Earth, two flybys of Venus, and three flybys of Mercury, a large propulsive maneuver 
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placed the spacecraft into a highly eccentric 12-h orbit around Mercury on 18 March 
2011. During the first year of orbital operations (primary mission), the spacecraft peri- 
apsis altitude was maintained between 200 km and 500 km with regular orbit-correction 
maneuvers (OCMs). Spacecraft angular momentum was controlled by internal reaction 
wheels and commanded momentum desaturations (CMDs). 

The periapsis latitude, initially at ~ 60°N, precessed northward to a maximum of 84.1°N 
in April 2013, after which the periapsis moved southward, reaching ~ 73°N by February 
2014. The maximum altitude, over southern polar latitudes, started at ~ 15, 000 km, 
a value that presents a challenge to the recovery of gravity anomalies in the southern 
hemisphere. After successful completion of the primary mission in March 2012 and in 
order to increase the frequency of observations at low altitudes, the spacecraft was placed 
in an 8-h orbit in April 2012, with an apoapsis altitude of ~ 10,000 km. With lower fuel 
reserves available, the periapsis altitude was allowed to evolve naturally and drifted to 
higher altitudes, reaching a maximum of ~ 450 km in early March 2013. Periapsis altitude 
is now decreasing progressively with each orbit, and an end of mission by impact onto 
Mercury’s surface is planned for March 2015. Figure 1 summarizes the orbit evolution of 
MESSENGER during the orbital phase of the mission. 

1.3. Outline 

In Section 2, we describe the available radiometric tracking data as well as the methods 
used to both process the data and obtain a gravity field solution. In Section 3, we introduce 
the new HgM005 gravity held and describe its geophysical attributes and implications. 
We then discuss in detail the associated sensitivity, uncertainties, correlations, and error 
calibration (Section 4). We also demonstrate how we used the MESSENGER range data 
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to improve the ephemeris of Mercury and the experimental determination of the 3 : 2 
resonance (Section 5). Finally, we show in Section 6 that near-circular frozen polar orbits 
are predicted to exist with the HgM005 gravity field. Such orbits could prove helpful to 
future exploration efforts because of lower orbit maintenance requirements and uniform 
global coverage. 

2. Data and Methods 
2.1. Data 

The MESSENGER spacecraft telecommunication subsystem operates in X-band (uplink 
at 7.2 GHz and downlink at 8.4 GHz) [Srinivasan et ai, 2007], with a typical noise 
level equivalent to 0.1 mm/s over a 60-s integration period. Two-way and three-way 
radio tracking data are acquired by the NASA Deep Space Network (DSN). The radio 
signals are quite sensitive to plasma noise, which affects the measurement noise level. 
The closer the radio signal path approaches the Sun (i.e., near superior conjunction, when 
Mercury passes behind the Sun as viewed from Earth), the more short-lived, turbulent 
plasma heterogeneities produce unknown shifts in the signal frequency received by ground 
stations. In those geometries, with small Sun-probe-Earth (SPE) angles (< 40°), the 
signal quality can be severely degraded to the point of becoming unusable for gravity field 
determination. 

The severe thermal environment at Mercury places important constraints on the oper- 
ation of the spacecraft. A fixed ceramic-cloth sunshade protects the spacecraft bus from 
solar radiation and limits possible spacecraft orientation by requiring the Sun direction 
to be within 10° of the sunshade-normal vector. Eor this reason, a variety of fixed an- 
tennas are used [Srinivasan et ai, 2007]. Each phased-array high-gain antenna (HGA) 
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is collocated with a fanbeam medium-gain antenna (MGA) to provide coverage in oppo- 
site hemispheres for different Earth positions relative to Mercury. There are four low-gain 
antennas (EGAs) that can provide coverage when the spacecraft performs specific observa- 
tions, and those are most commonly used near periapsis. For high-data-rate transmission 
passes, the radio-frequency system uses the fanbeam in uplink and the HGA in downlink. 
During these passes, the high radiometric signal levels reduce the receiver thermal noise, 
thus providing higher- quality tracking data (Doppler and range). Figure 2 shows that 
at favorable (large) SPE angles, the two HGAs have substantially better performance 
than the EGAs. At low SPE angles, the plasma noise dominates. Because of pointing 
constraints during operations, the HGAs were used mostly at high altitudes for downlink. 

Near periapsis, the spacecraft must maintain the main suite of instruments in a near- 
nadir orientation. During slews, an apparent line-of-sight velocity of the antenna is in- 
duced by the rotation of the spacecraft. Initially we found that Doppler residuals could 
show high-frequency patterns, just before and after closest approach. Further analyses 
allowed us to relate these patterns to rapid spacecraft slews, and to recognize a discrep- 
ancy in the position of the spacecraft center of mass (COM) reported by the spacecraft 
guidance, navigation, and control (GNC) team compared with the reference frame of the 
survey of the antenna phase centers. We used the Doppler data to identify this error of 
~ 90 cm in relative position and further refine the center-of-mass position. Figure SI 
shows the estimate of the center-of-mass position during the MESSENGER orbital phase. 
We adjusted the relative position of the center of mass and the antenna offsets after major 
OCMs, because of the larger expected shifts in COM position and absolute errors of the 
GNC reconstructions. Our solution is fully consistent with that of the GNC team, but it 
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is more accurate, with uncertainties at the mm level in each spacecraft reference frame 
direction. We have corrected the archived center-of-mass position information, and we 
note that the previous gravity solutions [Smith et ai, 2012; Mazarico et al, 2013] were 
not strongly affected by this error in the center-of-mass position, because we deleted the 
tracking data during the most rapid spacecraft slews. 

The gravity field presented here, HgM005, is an update to HgM002 [Smith et al, 2012], 
which was based on less than 6 months of orbital data, and to HgM004, a later solu- 
tion [Mazarico et al, 2013] that included an additional year of data. Here, we analyzed 
MESSENGER radiometric tracking data acquired by the DSN from the MESSENGER 
spacecraft through 4 February 2014. The minimum tracking altitude achieved in the or- 
bital phase versus position on Mercury is shown in Figure S2. In addition to measurements 
from the orbital mission phase, data from the three Mercury flybys are included. Notwith- 
standing the use of the tracking data from 2275 orbits spanning nearly 3 years, the data 
from the first two Mercury flybys still contribute to the determination of the equatorial 
gravity anomalies because of their low-altitude (~ 200 km) and near-equatorial closest 
approaches, compared with the ~ 1000 km equatorial altitudes during the orbital phase. 

A summary of the tracking data coverage and of the occurrences of spacecraft orbit 
maneuvers throughout the study period is given in Figure 3. MESSENGER was tracked 
by the DSN more extensively early in the mission. In the extended mission, only about 
one out of three periapsis passages was typically tracked, a schedule that yielded some 
one-day arcs with no data at altitudes sufficiently low to contribute to the gravity field 
determination. 
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For approximately two weeks around each of the seven superior conjunctions, the solar 
plasma noise was high. We did not include the weaker data from those arcs in the solution. 
This procedure was also followed for HgM002 [Smith et ai, 2012]. 

2.2. Methods 
2.2.1. GEODYN 

The orbit determination for MESSENGER has been performed using the GEODYN II 
orbit determination and geodetic parameter estimation software, developed and main- 
tained at NASA GSFC [Pavlis et al, 2013]. The MESSENGER tracking data were 
processed dynamically in 1-day segments (arcs), using a batch least-squares hlter [Mon- 
tenbruck and Gill, 2000; Tapley et al, 2004]. We explicitly modeled all forces acting 
on the spacecraft to integrate its trajectory, and we also modeled the radiometric data 
observables [Moyer, 2003]. The arc parameters were adjusted iteratively through least 
squares in order to yield the smallest observation residuals (i.e., discrepancies between ac- 
tual and predicted values). To limit the build-up of process noise and mismodeling errors, 
we reduced the arc length to no more than one day, i.e., two to three orbits, because the 
non-conservative forces such as that from direct solar radiation are large. MESSENGER 
experiences a solar flux that varies between 6.3 and 14.5 kW/m^ over the Mercury year, 
a variation due to the high eccentricity (~ 0.21) of Mercury’s orbit. For consistency, the 
arcs start and stop near apoapsis, which is convenient as it results in either two (before 
April 2012) or three MESSENGER orbits per arc. The full study period, March 2011 to 
February 2014, was divided into a total of 1058 one-day arcs. 
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2.2.2. General Models 


GEODYN relies on a number of models to integrate the spacecraft trajectory (force 
models), and to determine a computed observation to be compared with the actual ob- 
servable (measurement models). Here, we followed the approach outlined by Smith et al. 
[2012]. Briefly the measurement models include: troposphere refraction delays obtained 
from in situ meteorological data at the DSN sites. Earth orientation data supplied by the 
International Earth Rotation Service (lERS) [Gamhis, 2004], ocean loading corrections 
for the DSN sites obtained with the GOT4.7 ocean tide model [Ray, 2013, Appendix A], 
transformations between coordinate and atomic time [Moyer, 1981a, b], and spacecraft 
antenna offset corrections obtained from the spacecraft orientation and center-of-mass po- 
sition. The DSN observables have been modeled following Moyer [2003]. The force models 
include the gravitational acceleration associated with the gravity field of Mercury as given 
by the planetary gravitational constant, GM , the spherical harmonic coefficients, and 
Si^ for degree I and order m [Kaula, 1966], and the orientation model of Mercury; the 
relativistic modification of the central-body term; the third-body perturbations from ma- 
jor Solar System bodies computed from the DE423 planetary ephemeris [Folkner, 2010]; 
modeling of non-conservative forces such as those from solar radiation pressure and plan- 
etary radiation pressure using a spacecraft macromodel; and gravity tidal accelerations as 
predicted by the Love number k 2 - The modeling of planetary radiation pressure includes 
that due to Mercury’s albedo (reflected sunlight) and the planetary thermal radiation 

[Mazarico et al, 2012; Lemoine et al., 2013]. 

2.2.3. Non- conservative Force Models 
The non-conservative forces acting on the MESSENGER spacecraft must be modeled 

accurately in order to recover Mercury’s gravity field, tides, and rotational parameters. 
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The three major force models that limit the recovery of the pole orientation and tidal 
Love number in particular are solar, planetary albedo, and planetary thermal radiation 
pressures. As did Smith et al. [2012], we assumed a uniform surface albedo (0.074) and 
used a model of Mercury’s surface temperature [Paige et al., 2013] expanded to harmonic 
degree and order 4 to compute the surface thermal emission as a function of local solar 
time and latitude. However, because of its magnitude, the most important effect during 
the inversion results from the treatment of solar radiation. 

The radiation pressure model consists of a shape model with twelve plates: three plates 
represent the sunshade, hve plates the spacecraft bus (neglecting the panel behind the 
sunshade, as it is always occulted), and four panels for the front and back sides of the two 
solar panels. Telemetered quaternion data were used to orient the plates that represent 
the spacecraft bus, sunshade, and solar arrays. The solar panels can rotate independently 
from the bus, and we modeled their orientation around a gimbal as a function of time. 
This model explicitly includes the specular and diffusive reflectivity coefficients for each 
plate [Marshall and Luthcke, 1994]. 

Although we typically adjust a single parameter per arc to scale the three radiation 
accelerations, we find that the mismodeling of the radiation pressures can be better ac- 
commodated during arc convergence by estimating the areas of the 12 spacecraft plates. 
These areas generally increase the modeled spacecraft area by ~ 15-20% when the MES- 
SENGER orbit is in a noon-midnight configuration. In many instances, it is reduced, 
by up to ~10%, presumably due to self-shadowing effects not computed here [Mazarico 
et al., 2009]. During the global inversion of the solution, we kept the areas fixed to the 
values adjusted during orbit determination, but we estimated the sunshade and solar 
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panel reflectivities. The corrections to the a priori values are small, except for the diffuse 
reflectivity of the solar panel. We attribute this larger magnitude to the large thermal 
reradiation this panel must emit because of its high temperature, an effect not yet mod- 
eled in our work for MESSENGER, unlike the treatment by Antreasian and Rosborough 
[1992] and Marshall and Luthcke [1994] for TOPEX/Poseidon. The estimation of the 
plate areas during the orbit determination process and of the reflectivities in the global 
iteration help alleviate such model shortcomings. When spacecraft areas and reflectivities 
are not corrected, the radiometric data hts worsen, and the estimates of the Love number 
^2 and the obliquity are affected. Some of the changes in the low-degree held compared 
with HgM002 can also be attributed to the improved force modeling presented here. 

2.2.4. Planetary Orientation 

A gravity held model is inextricably linked to the dehnition of a model for planetary 
orientation. In the absence of a good orientation model, the position of gravity anomalies 
would not appear hxed over several Mercury rotations, and the anomaly estimate would 
be compromised. The orientation parameters can be considered force model parameters, 
as they affect the inertial trajectory of the satellite, which is the basis of estimation. The 
orientation parameters are also critical to ascertain the state of the interior of the planet 
(Section 3.6). Margot et al. [2007] determined an orientation model of Mercury from 
ground-based radar measurements, adding a non-zero obliquity and periodic longitudinal 
librations to the existing International Astronomical Union (lAU) model [Davies et al., 
1980]. Margot [2009] and Margot et al. [2012] further rehned those parameters, with a 
libration amplitude at the equator of ~ 450 m. We used their prime meridian to dehne 
the principal axes (PA) frame in which the gravity held is best expressed, rather than the 
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lAU prime meridian that was chosen to maintain the Hun Kal crater at a longitude of 
— 20°E. In this PA frame, the C 21 , S 21 , and ^22 coefficients are expected to be close to 
zero. 

For messenger’s orbital mission phase, the project adopted the JPL DE423 
ephemeris [Folkner, 2010] for mission planning and archived data products, and the grav- 
ity fields produced from MESSENGER data have been based on arcs reconstructed with 
that ephemeris. Recently, Folkner et al. [2014] produced a new Solar System ephemerides 
solution, DE430. Although we did not use the newer solution to produce HgM005, we 
evaluated it with the range data, and we found a clear improvement in terms of range 
residuals. In Section 5, we discuss the Mercury ephemeris and its estimation. 

2.2.5. Solution Strategy 

Each converged arc was processed in GEODYN to create normal equations, which were 
used to produce the gravity solution. In addition to arc-specific parameters, each arc’s 
normal equation contains the 2,597 Stokes coefficients {Ci^ for a degree and order 50 
gravity held, but also a suite of other global (common) parameters: the spacecraft 

low-gain antenna position correction, the spacecraft rehectivity parameters, the tidal Love 
number ^ 2 , the Mercury orientation parameters (right ascension, declination, and spin 
rate), and the Mercury gravitational parameter (GM). The antenna location adjustments 
are then tied together to correspond to a shift in the estimate of center-of-mass position. 

For the global solution, each converged arc was weighted according to its post- 
convergence observation residuals. In intermediary global solutions (predecessors to 
HgM005), arcs were weighted ~ 30% lower than their root mean squared (RMS) level; 
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the factor for HgM005 was 10%. This weighting scheme allows arcs of varying quality to 
contribute to the solution appropriately. 

The HgMOOS gravity field is based on 632 orbital arcs, out of the original 1058 arcs 
covering the study period. The former figure constitutes a substantial down-selection, 
but we found it beneficial to the overall solution. Although Smith et al. [2012] could not 
be as selective for HgM002 because of the short data span, now that we have about three 
years of orbital tracking data, we can select only the best-determined and most sensitive 
arcs without sacrificing broad spatial coverage. 

Because of solar plasma effects, we did not consider 337 arcs with low SPE angle (< 40°). 
We also excluded all the arcs with any kind of spacecraft maneuver, specifically 132 arcs 
over the three years of orbital observations. The total number of arcs removed was 426, 
accounting for the fact that some maneuvers occurred near superior solar conjunctions. 
Four of the remaining arcs were also deleted due to higher-than-normal residual RMS 
(> 1 mm/s). 

The normal equations were formed from each individual arc with GEODYN, and then 
merged with a degree power law (or Kaula) constraint equation before inversion using 
GEODYN’s companion program SOLVE [Pavlis et al, 2013]. The constraint is necessary 
to ensure the stability of the solution, given that the expansion degree is larger than what 
the data alone can support globally (especially in the southern hemisphere); this point is 
discussed in detail in Section 4.1. Although we experimented with a variety of strategies 
(e.g., the choice of parameters to estimate), the solution presented here uses an inversion 
approach for the global parameters similar to Smith et al. [2012]. Specular and diffuse 
reflectivity parameters for the solar arrays and specular reflectivity for the sunshade were 
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carried through to the final least-squares inversion; the other arc parameters were back- 


substituted during the combination of each arc’s normal equation. 

3. The HgM005 Solution, and Geophysical Implications 
3.1. Gravity Anomalies and Mascons 

The free air gravity anomaly held of the HgM005 gravity model is shown in map view 
in Figure 4, in both Mercator and polar stereographic projections. The map is visually 
very similar to those of Smith et al. [2012] and Mazarico et al. [2013] (Figure S3), despite 
the additional years of data. However, in HgM005 large-scale anomalies have greater 
conhdence, and additional shorter- wavelength signals are present because of the harmonic 
expansion to degree and order 50. 

The most prominent features in the gravity anomaly map (Figure 4) are the large 
positive anomalies, over the northern rise (~ 70° A, 35° E), over the Caloris basin (~ 30° A, 
160° A), and near the Sobkou basin (~ 35° A, 225° A). These features have been discussed 
by Smith et al. [2012]. As they argued, although only the Caloris anomaly is directly 
associated with an impact basin, the Budh-Sobkou anomaly may also qualify as a basin- 
associated mass concentration, or mascon, because it is a gravity high and is associated 
with an elevated crust-mantle boundary in crustal thickness models. Several apparently 
continuous linear features appear in Figure 4, such as one between (~ 10°A, 60°A) 
and (~ 30° A, 140° A). These linear features are correlated with topography, but they 
are relatively subdued in power, indicating that the topographic variations might be 
isostatically compensated. We present an updated crustal thickness map in Section 3.3. 
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3.2. Correlation with Topography 

As the Gravity Recovery and Interior Laboratory (GRAIL) mission demonstrated to 
extraordinarily high resolution at the Moon [Zuher et ai, 2013], gravity and topography 
are expected to be correlated, particularly at shorter wavelengths, because surface relief 
contributes to the planet’s gravitational potential. However, at the longest wavelengths 
(low degree and order), such a correlation does not necessarily hold because of isostatic 
compensation as well as uncompensated subsurface mass anomalies in the crust and man- 
tle. The gravitational signature of structures at depth is attenuated at orbital altitudes 
and can be detected only as long-wavelength signals. As on the Moon, in the spectral 
range of the HgM005 gravity held (harmonic degree 1 = 2 — 50), we do not expect the 
gravity and the topography to be fully correlated. Nonetheless, higher correlation values 
have typically been taken as indicators of gravity held improvement, in the case of Mars 
[Konopliv et al, 2011] and the Moon [Zuber et al, 2013]. 

We computed the global correlation of the gravity held with a global topography model. 
Smith [2014] archived in the NASA Planetary Data System (PDS) a spherical harmonic 
expansion of Mercury’s topography to degree and order 120 for a solution that combines 
Mercury Laser Altimeter (MLA) measurements [Zuber et al, 2012] and radio occultation 
measurements [Perry et al, 2013]. Because MLA cannot obtain measurements at ranges 
greater than ~ 1800 km, the occultation-derived radii, although limited in precision and 
spatial coverage, are critically important to dehne the long-wavelength shape of the south- 
ern hemisphere. However, in the northern hemisphere, the topography is not the limiting 
error source, as the density of MLA measurements is sufficient at the resolution of the 
gravity field in our analysis. As shown in Figure 5, the correlations for I < 20 (see Sec- 
tion 4.3) are typically ~ 0.6. This correlation is reasonably high given the challenges posed 
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by the eccentric orbit of MESSENGER. HgM005 has higher correlation values than the 
earlier HgM002 and HgM004 solutions, more markedly at / > 10 (Figure 5). For degrees 
greater than / = 20, the correlations steadily decrease to < 0.2 for /> 35, as expected 
from the degree strength (Section 4.3). This comparison indicates that the resolution of 
the gravity anomaly field in the southern hemisphere is low. This inference is further 
demonstrated in the lower panel of Figure 5, where we show localized correlations, com- 
puted following Wieczorek and Simons [2005]. These correlation values, after localization 
with a windowing width of either = 2 or l^in = 5 for a cap of half-angle 30°, are 
higher than the global values, typically 0.6 — 0.8 for / = 5 — 15. The higher correlation 
values over / = 25 — 35 are also indicative of better-resolved gravity anomalies for the 
northernmost latitudes. 

3.3. Bouguer Gravity Anomaly and Crustal Thickness 

With the same topographic model as that used in Section 3.2 for the correlation of 
gravity and topography, we computed the gravity expected from the topography. We 
assumed a uniform density for the crust p = 3200 kg m~^ [Smith et al, 2012], and we 
made use of the finite-amplitude method of Wieczorek and Phillips [1998] up to degree 
and order 5. 

The topography model has much higher intrinsic resolution than the gravity model, 
especially at lower latitudes. In order to be compatible with the spatial scales resolved 
by the HgM005 gravity field, we therefore limited the resolution of the topography field 
by expanding a spherical harmonic representation only up to low degrees. This procedure 
effectively filters out shorter-wavelength features. We found, however, that a single degree 
for truncating the expansion is not optimal. A truncation at / = 20 is adequate near the 
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equator, but such a limit is too severe at high northern latitudes, where HgM005 resolves 
shorter-wavelength structure. A truncation at / = 50 is more suitable for the north 
polar region, but it yields more features and power at lower latitudes than the gravity 
field can resolve. To construct the gravity field from topography shown in Figure 6a, 
we expanded the gravitational potential from topography to a degree consistent with 
the resolution of the gravity field at that location, as specified by the HgM005 degree 
strength (cf. Section 4.3). We then subtracted the topography-derived gravity from 
the measured HgM005 gravity held to obtain the HgMOOS Bouguer gravity anomaly held 
(Figure 6b). At the northernmost latitudes, several features in the free air gravity anomaly 
and Bouguer correction maps are not seen in the Bouguer anomaly held, indicating a lack 
of compensation. This effect is notably the case for the northern rise mentioned above. 

Under the assumption of Airy isostasy, the Bouguer anomalies are indicative of varia- 
tions in the depth of the crust-mantle boundary, and the Bouguer map can be translated 
into a map of crustal thickness (Figure 6c). Following Smith et al. [2012], we assumed an 
average crustal thickness of 50 km and a crust-mantle density contrast of 250 kg m“^. The 
equatorial and polar regions are characterized by thicker and thinner crust than average, 
respectively. Such long-wavelength variations in apparent crustal thickness might also 
have contributions from variations in crustal or mantle density. 

3.4. Low-degree Field 

As noted above, the low-degree coefficients in the gravity held are important both 
for understanding the structure of Mercury’s interior and for modeling the long-term 
evolution of a spacecraft orbit. The HgM005 C 20 and C 22 values are in good agreement 
with the HgM002 values of Smith et al. [2012], having changed respectively by 0.22% and 
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0.67%, figures less than the standard deviation of 0.01 x 10“® (i.e., 0.44% and 0.80%, 
respectively) in the HgM002 solution. The differences with respect to HgM004 are even 
smaller (0.05% and 0.27%, respectively), indicating the robustness of these estimates with 
increased temporal coverage. However, larger discrepancies (6 — 9%) exist for the C30, 
C40, and higher-degree zonal coefficients, which are highly correlated (Section 4.5). The 
HgM005 values are closer to the estimates of Genova et al. [2013] than to the HgM002 
values. We ascribe those changes to a more careful consideration of the antenna phase 
offset corrections and of the parameters associated with solar radiation (area scale factors 
and panel reflectivities), as described in Section 2.2.1. These changes are particularly 
relevant to the orbit evolution of the Mercury Planetary Orbiter of the BepiColombo 
mission now in development by the European Space Agency and the Japan Aerospace 
Exploration Agency, as discussed by Genova et al. [2013]. 

The gravitational parameter of Mercury {GM) also differs from previous estimates, at 
2.2031870799 x 10^^ ± 8.6 x 10® m® s“^ (after calibration, see Section 4.4). Although the 
changes from HgM002 and HgM004 (~ 9 x 10^ and ~ 3 x 10^ m® s“^ respectively) are small 
in absolute terms, they are significant compared with the formal uncertainties. Again, this 
improvement is likely due to improved (reduced) correlations with the solar radiation and 
the spacecraft state-vector parameters. The HgM005 value of GM is in good agreement 
with the navigation team determination (difference of ~ 1 x 10^ m®s“^), especially when 
compared with the estimates from Anderson et al. [1987] and Genova et al. [2013] (with 
differences of approximately —22 x 10^ and —19 x 10^ m® s“^, respectively). 
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3.5. Tidal Love Number ^2 

The tidal Love number ^2 describes the amplitude of the time- variable degree-2 gravity 
signal due to the tides raised by the Sun on Mercury. In addition to its direct effect on 
the trajectory of MESSENGER, it can provide constraints on the internal structure of 
Mercury [Van Hoolst and Jacobs, 2003; Van Hoolst et ai, 2007; Rivoldini et ai, 2009; 
Padovan et ai, 2014], beyond those imposed by the moment of inertia and obliquity 
[Peale et ai, 2002; Margot et ai, 2012; Ptauck et al, 2013]. The parameter is particularly 
sensitive to the core size and outer solid shell thickness, and it also varies with the rigidity 
and temperature of the mantle layer, as illustrated by Padovan et al. [2014]. The HgM005 
solution for the tidal Love number, k 2 = 0.451 ±0.014 (after scaling by a calibration factor 
of 10, see Section 4.4), is the first obtained directly from observations. We note that our 
value depends on the specific modeling and inversion strategy of the radiation pressure 
accelerations. Accounting for possible systematic effects and biases, a wider range of 
k 2 = 0.43-0.50 can thus not be ruled out with our HgM005 solution. 

We used the ALMA modeling tool [Spada, 2008] to calculate the k 2 value expected from 
models of internal structure. Following Smith et al. [2012] and Hauck et al. [2013], we 
varied the thickness of the lithosphere between 70 and 90 km and the rigidity of both the 
lithosphere and a possible solid FeS layer (between the mantle and the fluid core). The 
resulting k 2 values (Table 1) range between 0.46 and 0.62. The range is consistent with 
that of Van Hoolst and Jacobs [2003] and Rivoldini et al. [2009], who suggested a value 
of 0.4-0. 6 prior to the MESSENGER mission. We note that the larger values (/c 2 > 0.5) 
are obtained only when assuming that a solid FeS layer is present. Initially proposed by 
Smith et al. [2012] to account for the large moment of inertia of the solid outer shell, as 
noted above, an updated analysis by Hauck et al. [2013] with the latest obliquity values 
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[Margot et al, 2012] found that such a layer, although still compatible with the data, is 
no longer required, in agreement with the findings of Rivoldini and Van Hoolst [2013]. 

More recently, Padovan et al. [2014] explored the tidal Love number with a rheological 
model and a large range of governing parameters to hnd solutions compatible with the 
MESSENGER and ground-based geophysical observations [Hauck et ai, 2013]. Although 
somewhat extreme scenarios with very high rigidities can produce low k 2 values as low as 
0.4, Padovan et al. [2014] obtained a range of 0.45-0.52 for a mantle grain size of 1 cm. 
Smaller grain sizes (1 mm) would result in larger values (by ~ 10%). Our estimate of /c 2 
from the gravity field solution is fully consistent with their results, leaning slightly toward 
the lower values. 

The k 2 value obtained from MESSENGER tracking data, although still considered pre- 
liminary, is fully consistent with model expectations built from the earlier gravity results 
of Smith et al. [2012]. As additional gravitational data are acquired by MESSENGER, 
increasingly rehned estimates of k 2 will help to further constrain the interior structure 
and rheology of Mercury’s interior. 

3.6. Mercury Orientation: Pole Position, Obliquity, and Spin 

With less than 6 months of data, i.e., less than 3 Mercury spin periods, Smith et al. [2012] 
did not attempt to estimate Mercury’s orientation parameters along with the HgM002 
gravity field. For the HgM005 solution, in contrast, we estimated the right ascension 
(RA) and declination (DEG) of the pole and the spin rate. We prepared an alternate 
solution for which we also adjusted the amplitude of the longitudinal librations, but our 
current sensitivity to that parameter with current radio tracking data is limited. 
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Because the pole parameters are connected to the C^i-, and S 22 terms, small es- 
timated values for these coefficients are indicative of good recovery of pole parameters, 
because the expected values of these coefficients in the PA frame are zero. As shown 
in Table 3, we do obtain suitably low values for those coefficients. We note that the 
C 22 and S 22 values may indicate a shift of the principal axes frame in longitude of 
5<t>PA = arctan{S 22 / C 22 ) /‘^ ~ —0.048°, equivalent to ~ 2 km at the equator. Use of 
the lAU convention for the prime meridian [Archinal et ai, 2011] would result in a sub- 
stantially larger value of S 22 - 

Our updated pole position at J2000 is RA = 281.00480° ± 0.0054 and DEC = 
61.41436° ± 0.0021 (after calibration, see Section 4.4), approximately 10 arcseconds away 
from the best-fit position of Margot et al. [2012] which we used as an a priori estimate in 
our solution. From these values, we computed the obliquity, the angle between Mercury’s 
orbit plane normal and its spin axis. Margot [2009] obtained an obliquity of 2.11 ± 0.10 
arcmin. With improved modeling and additional ground-based radar observations, Mar- 
got et al. [2012] refined the spin axis orientation {RA = 281.0103°, DEC = 61.4155°) and 
revised their obliquity estimate to 2.04 ± 0.08 arcmin. The HgM005 solution yields an 
obliquity of 2.06 ± 0.16 arcmin (calibrated uncertainty). This value is entirely consistent 
with both the Margot [2009] and Margot et al. [2012] estimates (within half a standard 
deviation) . 

From the C 20 and C 22 values and this obliquity (with an uncertainty equal to the 
formal error multiplied by 5), we obtain an estimate of the polar moment of inertia 
C/MR'^ = 0.349 ± 0.014, where R is Mercury’s mean radius. The fractional part due to 
the solid outer shell is CmjC = 0.424 ± 0.024. These values are only slightly different from 
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those of Margot et al. [2012] (larger by 0.2 and 0.25 standard deviations, respectively). 
Such small changes are compatible with the majority of the interior structure models 
presented by Hauck et al. [2013]. Compared with their nominal model, our updated 
values imply a slight decrease in the outer shell thickness (by ~ 15 km), an increase in 
the outer shell density (by ~ 54 kg m“^), and an increase in the density of the innermost 
solid and outer liquid cores (by ~ 115 kg m“^). 

With Mercury in a 3: 2 spin-orbit resonance, its spin rate (except for the effect of 
librations) is directly tied to Mercury’s orbital period (around the Sun). We obtain a spin 
rate correction relative to the lAU convention (a spin period of 58.646220 days) [Archinal 
et al, 2011] of (9.042 ± 1.288) x 10“^^ degree s“^ or a corrected period of 58.646146 ± 
0.000011 days. In addition to the good sensitivity of the MESSENGER observations over 
nearly 20 spin periods, this change is justihed by the improved consistency with the 3; 2 
resonance when combined with the orbital period (Section 5). The lAU convention does 
not include an uncertainty, but our new estimate is consistent with the spin period of 
Klaasen [1976] (58.6461 ± 0.005 days). In Section 5, we discuss the adjustment of the 
Mercury ephemeris, which can be used as a further check on this updated rotation period, 
because of the 3: 2 spin-orbit resonance. 

4. Sensitivity Analysis 

In this section, we evaluate the quality, sensitivity, and uncertainties of the HgM005 
gravity field. We first justify our choice of the regularization constraint to stabilize the 
gravitational anomalies in the southern hemisphere (Section 4.1). We then propagate 
the covariance matrix obtained during the least-squares inversion to spatially map the 
expected error levels in the anomaly map (Section 4.2), and we construct a degree strength 
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map describing the spatial resolution of HgM005 (Section 4.3). We determine a scale 
factor to calibrate the formal uncertainties to more conservative values (Section 4.4), and 
we evaluate the correlation between estimated parameters, in particular the zonal terms 
(Section 4.5). 

4.1. Kaula Constraint 

Regularization is required since the solution is being expanded to a spherical harmonic 
degree that the data do not support on a global basis, because of a strong variation of 
altitude with latitude. 

An a priori constraint on the magnitude of the Stokes coefficients is necessary for the 
determination of spherical harmonic gravity field solutions if the data are not globally 
distributed at the wavelength of the truncation degree. In the determination of planetary 
gravity fields, the “Kaula rule” is used for a smoothing constraint, whereby each coefficient 
or §1^ is assigned an a priori uncertainty on the basis of its expected variance at 
degree I [Kaula, 1966]. Empirically, this constraint follows a 1//^ relationship, which can 
be justified a priori by the self-similar fractal nature of planetary surfaces and a posteriori 
by the gravity power spectra of Earth [Lemoine et ai, 1998], Mars [Lemoine et al, 1997], 
and the Moon [Lemoine et al, 2013]. The strength of the constraint is given by a scaling 
factor K, for an expected RMS power of K/t^. Each planetary body has its own Kaula 
scale factor. Through scaling relationships from other bodies [Konopliv et al, 2014], we 
would expect K 5 X 10 ® for Mercury (Table 2). 


ure 7 does match the TC ~ 4 x 10 ® rule reasonably well at low degrees, indicating that 
the “true” level of the Kaula factor may be close to that number. Unfortunately, because 
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The RMS power 






of messenger’s eccentric orbit, use of this Kaula factor globally leads to implausible 
gravity anomalies in the southern hemisphere. We therefore chose a stronger Kaula con- 
straint, K = 1.25 X 10“®, as it prevents large anomalies from developing in the southern 
latitudes without overly smoothing the solution in the north. We found that strength- 
ening the constraint further is detrimental, as it smoothes the field nearly equally in the 
south and in the north. 

Because a consequence of the Kaula constraint is to bias each coefficient toward zero, we 
prepared low-degree solutions with no Kaula constraint in order to ascertain that the low- 
degree field is generally not affected by its application. The low-degree field is especially 
important for interior modeling [Hauck et ai, 2013] and calculations of spacecraft orbit 
evolution [Genova et al, 2013]. We obtained alternate solutions with no Kaula constraint, 
and with increasing truncation degrees up to / = 20. These unconstrained solutions are 
themselves strongly degraded, because of truncation aliasing (for the I = 4-5 solutions) 
and instabilities leading to large power at the higher degrees (for the I = 7-20 solutions). 
Nonetheless, the RMS power of these unconstrained solutions indicates that the HgM005 
solution is not driven to a lower power because of the Kaula constraint (Figure 7). As an 
additional precaution, we applied the Kaula rule only above I = 7 for HgM005, whereas 
Smith et al. [2012] applied the Kaula constraint starting at degree I = 3. However, we 
note that the differences from a solution constrained from I = 2 are very small, with the 
RMS power of the differences lying below the formal error spectrum. These unconstrained 
solutions also show that without regularization the MESSENGER data can determine a 
gravity field only to degree and order 6. Gomparison of gravity anomalies expanded to 
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I = 15 with those obtained with a stronger constraint {K = 4.0 x 10“®) shows small 
differences in the northern hemisphere, below 2 mGal RMS. 

At high degrees (/ > 20), the power drops markedly, to approach the level of the 
K = 1.25 X 10“® rule. We interpret this behavior to signify that the data can partially 
support a degree I = 15-20 field, as the power up to / = 20 stands high despite the chosen 
Kaula rule. The Kaula constraint drives the global power levels at shorter wavelength 
(/ > 20), because the data are sensitive only to anomalies at the northernmost latitudes. 
Thus, we select K = 1.25 x 10“® for its effect at high degrees, where it is necessary, rather 
than for its validity at long wavelengths, where it is less needed. More refined approaches 
exist and could be explored for Mercury, such as spatial-spectral constraints [Konopliv 
et al, 1999] or localized Kaula constraints [Han et ai, 2009; Mazarico et ai, 2010], but 
such efforts are beyond the scope of this paper. 

As we explained above, we recognize that the total field power may be over-constrained 
by the stronger Kaula factor, and we expect the true gravity field power to follow a 
A" ~ 4 X 10“® rule. We justify our choice because for the primary purposes of HgM005, 
including global geophysical analysis, we want to prevent large but poorly constrained 
anomalies in the southern hemisphere. Future low-altitude data from MESSENGER or 
other spacecraft will be important in refining the gravity field. 

4.2. Covariance Error Analysis 

The projected formal errors for the gravity anomaly field for the HgM005 model, as 
calculated from the full error covariance matrix, are shown in Figure 8a. The formal 
errors range between 2.5 and 8.5 mGal and vary predominantly with latitude because of 
the eccentric orbit of MESSENGER. Comparison with earlier helds (HgM002, HgM004) 
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shows typical 2-4 mGal RMS differences northward of 20°N, in good agreement with the 
propagated errors. In the southern hemisphere, RMS differences range from 6 to 18 mGal 
with an average of 10 mGal, indicating that the formal errors underestimate uncertainties 
at those latitudes. The low-altitude flyby passes positively influence the projected errors 
near the equator. We also calculated the errors from only those coefficients up to / = 20. 
The northern hemisphere anomaly errors at those longer wavelengths are much smaller, 
a result that is consistent with the error power spectrum in Figure 7 reaching a noise- 
to-signal ratio near unity around / = 25 — 35. In the south, where the Kaula constraint 
is more effective in reducing the power of the shorter wavelengths, the error level is not 
decreased as appreciably and is still much higher than in the north. 

Taking into account the calibration factor discussed below, it might be possible that 
gravity anomalies as large as 60-80 mGal, and not confined to short wavelengths, could 
remain undetected at southern latitudes. 

4.3. Degree Strength 

An extension of this spatial covariance error propagation is the degree strength map, 
illustrated for Venus by Konopliv et al. [1999] and for the Moon by Konopliv et al. [2013]. 
The idea is to obtain, for each point on the globe, the degree at which the anomaly error 
calculated from the covariance matrix is equal to the expected anomaly given that the 
gravity coefficients follow the Kaula constraint precisely, that is, when the signal-to-noise 
ratio at that point is equal to unity (as illustrated in Figure S5). Clearly, the resulting 
map, shown in Figure 8b, has a zonal pattern because of MESSENGER’S eccentric orbit, 
much like the anomaly error map. The minimum degree strength is found in the south 
polar region, where htrength ~ 8, and the maximum Istrength ~ 36 is near the north pole. 
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Near the equator, the degree strength is ~ 15 corresponding to a resolution of ~ 500 
km on the surface, which is perhaps optimistic despite the two short tracked segments of 
the flybys given that the spacecraft altitude over the equator is generally ~ 1000 km. A 
degree strength of 35-36 near the north pole is in line with the expected resolution given 
spacecraft altitudes above 200 km. 

The degree strength is useful in gauging the extent to which the anomaly map can be 
used with conhdence for geophysical interpretation at a given location. To illustrate this 
point, we used the degree strength map to create an anomaly map expanded at each point 
only up to its specihc degree strength. Figure 8c presents this “degree strength anomaly 
map,” which is smoother than the full expansion of HgM005 (Figure 4) in the northern 
hemisphere and slightly damped near the equator. However, it is better resolved than the 
I = 20 expansion of HgM005 (and HgM002), as its degree strength reaches ~ 36. 

4.4. Error Calibration 

Smith et al. [2012] scaled the formal uncertainties for HgM002 on the basis of their 
understanding of the confidence level that could be assigned to the low-degree coefficients. 
In particular, the quoted 0.01 x 10“® uncertainty for the I = 2 coefficients reconciled the 
magnitude of the C 21 , S 21 , and S 22 coefficients with their expected zero values in the 
principal axes frame. We use the same argument (Table 3). 

The HgM005 low-degree formal uncertainties are smaller by a factor of ~ 5 than those for 
HgM002. We also consider the magnitude of the gravity anomaly error implied spatially 
by the HgM005 covariance matrix (Section 4.2): given the existence of ~ 120 mGal 
gravity anomalies in the northern hemisphere, and the expectation that some unresolved 
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or indiscernible anomalies of similar magnitude might exist in the southern hemisphere, 
uncertainties of 80 niGal or more in that region are appropriate. 

As a result, we recommend a scale factor of 10 to 15 to obtain conservative error 
estimates for the gravity field coefficients. This recommendation translates to a scaled 
uncertainty of 0.001 x 10“® on the C 20 and C 22 coefficients, still a ten-fold improvement 
relative to HgM002. The tidal and rotational parameters, because of their long wavelength 
and global scale, do not suffer from the same loss of sensitivity because of orbit geometry, 
and we do not apply a scaling of the formal uncertainties. 

4.5. Low-degree Coefficient Correlations 

As noted by Smith et al. [2012], the low-degree zonal coefficients are highly correlated. 
Due to the nature of their perturbations, neighboring zonals are anti-correlated, but the 
correlation coefficients rapidly decrease with a mismatch between degree and/or order 
(see Figure S4). The slow rotation of Mercury and the short arc duration chosen because 
of force mismodeling concerns make the sampling of the resonances due to the zonal 
terms difficult (for example, the spacecraft initial state is estimated daily), and we can 
observe only a lumped effect. This effect is, however, not limited to zonals, as or Si^ 
coefficients of the same order m exhibit this behavior, which is expected as such terms 
produce perturbations at the same frequencies [Kaula, 1966]. 

The additional data included in HgM005 compared with HgM002 did not alleviate this 
behavior, but the improved force modeling and inversion strategies of arc selection and arc 
weighting helped reduce the strongest correlations. For example, the correlation between 
C20 and C30 decreased from —0.86 to —0.67. This decreased correlation adds confidence 
to the determination of the C30 coefficient, important for long-term orbit evolution (Sec- 
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tion 6). Future improvements may be obtained when performing the orbit determination 
over longer arcs, and of course once other spacecraft are placed into different orbits around 
Mercury, such as BepiColombo [less et ai, 2009]. 

5. Mercury’s Ephemeris and 3: 2 Resonance 

The range measurements to MESSENGER, although not contributing as much as the 
range rate to the gravity held solution, are directly sensitive to errors in the modeled 
position of the planet, or ephemeris error. Figure 9 shows the pass-by-pass range biases 
estimated after convergence of the one-day arcs with the MESSENGER radio tracking 
data given the starting DE423 ephemeris. The periodic variations indicate Mercury posi- 
tion errors. Whereas the DE423 and DE430 planetary ephemerides result from a combined 
trajectory adjustment of the full Solar System (i.e., hundreds of major and minor plan- 
etary bodies), in this work we attempt only to adjust the orbital elements of Mercury 
itself. With the more recent DE430 ephemeris (which includes some MESSENGER data) 
[Folkner et al, 2014]), the higher-frequency errors are signihcantly reduced: the RMS of 
the range residuals outside of the low-SPE-angle periods decreases from 57.8 m to 24.4 m 
(Figure 9). 

After having obtained the HgM005 solution, we re-converged the data arcs with our new 
gravity held. We used the residuals of the range data over the full mission to estimate a 
relative correction at J2000 to the trajectory of Mercury, with the “Set III” formulation 
of Brouwer and Clemence [1961] (Table 4). We then reprocessed the data arcs with 
these corrections. We hnd signihcant improvements with this linear correction in J2000 
compared with both the starting ephemeris DE423 and the DE430 solution, with a range 
residual RMS of 10.9 m (Figure 9). In-plane corrections are ~ 50 m. Out-of-plane 
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corrections are larger, about 1 km, but we note that the cross-track (inclination) direction 
is not as well constrained from ground-based radio tracking. The incorporation of a long 
time series of Earth-MESSENGER range measurements will likely improve the quality of 
the Solar System and Mercury ephemerides further. We note that Verma et al. [2014] 
recently improved their Solar System ephemeris using MESSENGER tracking data and 
obtained a similar level of post-fit range residuals (2.8 ± 12.0 m with data up to September 
2012 ). 

Our correction also yields an updated semi-major axis of Mercury, which can be con- 
verted to an orbital period around the Sun, to which we add ~ 8 s to account for the preces- 
sion of the longitude of the perihelion [Shapiro, 1989]. We hnd Porbu = 87.969216879 days 
± 6 s. In combination with our new spin period estimate of Pspin = 58.646146 ± 0.000011 
days (Section 3.6), we find a ratio PorUt/ Pspin of 1.49999900. This ratio is a factor of ~ 2 
closer to the expected 3: 2 ratio than for our a priori models (the orientation parameters 
of Margot et al. [2012] and the DE423 ephemeris), with the discrepancy decreasing from 
~ 3.5 X 10“® to ~ 1.0 X 10“®, and it provides further and more accurate experimental 
determination of Mercury’s lock in the 3: 2 resonance. It also demonstrates the quality of 
our two independent measurements of Mercury’s spin rate and ephemeris. 

6. Future Mission Planning, and BepiColombo 

With confidence in the low-degree gravity field, it is possible to consider whether there 
are certain classes of orbits that are particularly advantageous for future Mercury orbital 
missions. Although planning for the next mission to Mercury, the dual BepiColombo 
orbiters, is too far advanced to beneht from this analysis, the exercise is more than of aca- 
demic interest. Indeed, the MESSENGER mission lifetime was restricted by the amount 
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of propellant available after orbit insertion. So-called “near-frozen” orbits, defined as or- 
bits for which the elements are nearly constant when averaged over a finite time interval, 
offer a cost-saving alternative to the frequent orbit-correction maneuvers MESSENGER 
employed to manage periapsis altitude and delay impact. 

Previous work on frozen orbits around Mercury include the analyses by Debate et al. 
[2010] and Ma and Li [2013], but we note that those studies did not benefit from or make 
use of the recent MESSENGER gravity results. They also considered only the C 20 and 
C 30 coefficients, although we note that they focused on highly eccentric orbits. In those 
cases, the perturbations by the gravitational attraction of the Sun dominate, along with 
the secular changes resulting from C 20 and C 30 . Here, we concentrate on near-circular 
frozen orbits, following Cook [1991], and we initially consider all the zonal coefficients 
of HgM005, i.e., up to degree 50. We perform a search for frozen orbits for different 
semi-major axis values over all values of inclination i from 0° to 180°. For each case, we 
obtain a frozen eccentricity. We find that the frozen eccentricity converges after inclusion 
of terms through degree I = 15 — 25, depending on the altitude considered. This relatively 
high-convergence degree may be a consequence of the high correlation of the zonal terms 
(Section 4.5). We note that the secular effects of the zonal terms decrease with degree 
and altitude, and that the power of the HgM005 zonals at high degree is dampened by the 
Kaula constraint. Nonetheless, it appears that C 20 and C 30 are not sufficient to predict 
the frozen eccentricity. 

The results of our calculations for an average altitude h of 300, 500, and 1000 km are 
show in Figure 10. The maximum allowable eccentricity (i.e., the eccentricity that would 
lead to a periapsis equal to the reference radius R = 2440 km) is, respectively, Cmax = 0.11, 


@2014 American Geophysical Union. All Rights Reserved. 



0.17, and 0.29. Frozen orbits exist for many of these inclinations and semi-major axis 
values. Outside of a narrow range (i = 65 — 70°), the periapses are typically near the 
south pole (argument of pericenter u = 270°), as in the case of the Moon (for instance, 
the orbit chosen for the Lunar Reconnaissance Orbiter during its commissioning and 
extended mission phases [Chin et ai, 2007]). Over a wide range of moderate inclinations 
(z = 30 — 60° prograde and retrograde), low-eccentricity frozen orbits are possible. For 
h = 1000 km, near-equatorial and near-circular frozen orbits exist. However, we focus 
here on the polar orbits, which are most interesting for a potential future orbiter mission, 
as they would provide global coverage. 

Because the Sun exerts a large third-body perturbation at Mercury and dominates 
orbit evolution for MESSENGER, we used GEODYN to perform a high-fidelity orbit 
propagation of the h = 500 km polar orbit, with the full (degree and order 50) HgM005 
gravity field (in contrast to the initial search discussed above for which only zonal terms 
were considered). In addition to the gravitational acceleration from the Sun, we also 
considered additional perturbations such as the direct solar radiation pressure and the 
planetary thermal and albedo radiation pressure accelerations. Eigure 11a shows the 
evolution of this h = 500 km polar orbit over 8.5 years in terms of equinoctial elements, a 
phase space used to evaluate orbit stability: esin{ijj) versus ecos{uj). The stability of the 
orbital elements is clear despite their complicated pattern, and no long-term drift exists. 
The periapsis altitude is also very stable (Eigure 11b), varying only by ~ 12 km over the 
8.5-year orbit integration. 

In order to consider the impact of the HgM005 uncertainties in our calculations, we 
performed those computations again, with 25 so-called “clone fields.” Each clone repre- 
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sents an alternative solution to the least-squares inversion, consistent with the covariance 
matrix. Varying each coefficient according to its individual uncertainty would not be 
appropriate, as doing so ignores the correlations between coefficients and would thus not 
statistically replicate the initial held. Instead, we construct the clone helds from the 
covariance matrix, which precisely accounts for the error characteristics of the HgM005 
gravity model. As described in the supplementary material of Smith et al. [2012], the 
covariance matrix was hrst diagonalized. Then each eigenvector was scaled by the square 
root of its eigenvalue and a random factor, and hnally the eigenvector was added to the 
baseline solution. Smith et al. [2012] used a Rademacher distribution for the random 
factors of their 50,000 clones (only values of -|-1 or —1 were allowed). Here, we use a 
Gaussian distribution, which enables the smaller number of clones to better represent the 
range of variability. 

We perform the frozen eccentricity search at different inclinations with 25 clones of 
HgM005. The variability in frozen eccentricity reduces the suitable regions, in particular 
near i = 0° ,i = 65°, and i = 115° (see also Figure S6). Whereas the near-polar inclinations 
show more susceptibility to the HgM005 uncertainties than i = 15 — 45° (and retrograde), 
they show that polar frozen orbits exist for eccentricities near 0.07, the value found with 
HgM005. 

We performed propagations with GEODYN of polar orbits with eccentricities of 0.064, 
0.068 and 0.074 (the spread over the 25 clones for h = 500 km and i = 90°). Each 
orbit conhguration was integrated with HgM005 and with the two clone helds bounding 
the eccentricity values. Although the initial periapsis altitudes differ, of course, the or- 
bit evolution is rather slow and appears sufficiently robust for future mission planning 
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consideration (Figure S6b). The long integration times (8.5 years) show that the initial 
larger variations are themselves periodic (~ 6 years) and would not lead to an impact. 
Near-circular polar frozen orbits are thus likely to exist at Mercury, a hnding that would 
benefit the prospects for long-term exploration and monitoring of the inner planet and its 
dynamic environment. Of course, the orbit dynamics alone do not dictate the mission de- 
sign, and the thermal environment in particular could be the most important constraint. 
We compared the planetary thermal radiation acceleration received along the orbit shown 
on Figure 11 with that received by MESSENGER during its hrst 100 days in orbit. Al- 
though the maximum flux is not significantly larger than for MESSENGER (Eigure S7), 
the shorter orbital period (~ 2 h versus ~ 12 h) sharply reduces the ability of the space- 
craft to cool on the nightside. A more detailed study would be valuable but is outside the 
scope of this work. 

The implications of the improved value of the C30 coefficient are especially important 
in the case of the BepiGolombo Mercury Planetary Orbiter, as the initial orbit design 
anticipated the (then unknown) C30 to be positive. Its estimated value, (—0.47659 ± 
0.0016) X 10“® (Table 3), thus leads to a decrease of the periapsis, and substantially 
more thermal forcing on the spacecraft components. This heating could be remedied by 
increasing the initial orbit altitude, and the better determination of C30 will facilitate the 
necessary modihcations in mission planning. 

7. Summary and Conclusions 

We have analyzed three years of radio tracking data collected at Mercury by the MES- 
SENGER spacecraft. We obtained a gravity field solution expanded to spherical harmonic 
degree and order 50, called HgM005, for which we also estimated the planetary orientation 
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and ephemeris of Mercury. We have described the geophysical implications of these new 
results, and we discussed in detail the modeling and error sources associated with the 
gravity anomalies and other important gravity parameters such as the low-degree zonal 
harmonics. 

After successfully completing its one-year primary and one-year hrst extended missions, 
sufficient fuel reserves remained on MESSENGER to design a novel end-of-mission sce- 
nario. MESSENGERS second extended mission will take advantage of the decrease in 
periapsis altitude due to solar perturbations to execute four low-altitude campaigns, each 
spanning several weeks with periapsis altitudes lower than 100 km and as low as 25 km. 
During the hrst two of these periods, in September and October 2014, the periapsis will 
be in view of the Earth. With a periapsis latitude as low as 65° the tracking of the 
spacecraft by the NASA DSN will give an exceptional view of the short-wavelength grav- 
ity anomalies over a large part of the northern hemisphere. For example, the majority 
of the western hemisphere between 50° A^ and 75° N will be mapped from altitudes less 
than 50 km. These data will yield a degree strength above 50 in the northern hemisphere, 
with the potential to substantially improve our understanding of the crustal structure of 
Mercury. 
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Figure 1. MESSENGER periapsis altitude (blue), periapsis latitude (red), and Sun-probe- 
Earth (SPE) angle (green) during the orbital mission phase, including predictions for the remain- 
der of the second extended mission (XM2). 
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Figure 2. Noise level of the individual tracking passes of MESSENGER, measured as the root 
mean square of the Doppler residuals over short and detrended segments and plotted against 
the Sun-probe-Earth angle. Each integer on the color bar indicates a separate MESSENGER 
antenna, following the PDS Frames Kernel for the radio science experiment. Antennas 0 and 3 
are HGA/MGAs, whereas the rest are EGAs. 
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Figure 3. Summary of tracking data coverage of the MESSENGER spacecraft during the 
orbital mission phase. Each line represents one orbit, and DSN passes are shown in red. Black 
dots indicate a spacecraft maneuver. 
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Figure 5. Correlation of MESSENGER-derived gravity models with the global shape model 

determined from ML A and radio occultation measurements. We show the global correlations 

(top) and the correlations of the fields after localization around the north pole with different 

windowing tapers (bottom). K is the Kaula factor, with a smaller number indicating a stronger 

constraint on the field power (see Section 4.1). 
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Figure 6. (top) Gravity anomaly field (in mGal) predicted by surface topographic relief, in- 
ferred from MLA measurements and radio occultations, for a crustal density of p = 3200%m“^. 


At left is a Mercator projection to 67° latitude, and at right a north polar stereographic pro- 
jection for latitudes 60-90°N. (middle) Bouguer gravity anomaly held (in mGal), obtained by 
subtracting the top gravity held from the free-air gravity anomaly held in Figure 4. (bottom) 
Crustal thickness (in km), obtained from hrst-order downward continuation of the Bouguer grav- 
ity anomaly. At each location, the spherical harmonic representation is expanded only up to the 

spatial degree strength of HgM005 (Figure 8). 
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Figure 7. RMS power of gravity solutions obtained with different Kaula factors (thick colored 
lines) and associated formal error spectra (thin colored lines). Our preferred solution, HgM005, 
is shown in red. The a priori Kaula factor {K = 4 x 10“®) is shown by a thin dashed line. The 
preferred Kaula factor {K = 1.25 x 10“®), shown as a thick dashed black line, was applied for 
I > 7. 
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Figure 8. (top) Gravity anomaly errors (in mGal) obtained from the HgM005 covariance, 
(middle) HgM005 degree strength computed from the comparison of the Kaula regularization 
and the full covariance matrix. Features smaller than the degree-strength-equivalent wavelength 
are not robust, (bottom) Free-air gravity anomaly field (in mGal) of HgM005, expanded at every 
location only up to the local degree strength. Same projections as in Figure 6. 
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• DE423: 0.87 ± 56.69 m 

• DE430: 2.12 ±24.12 m 

• this work: -0.16 ± 11.01 m 
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Figure 9. Range residuals during MESSENGER’S orbital mission phase, obtained after arc 
convergence with the DE423 ephemeris (green), the DE430 ephemeris (red), and our adjustment 
(blue). Each point represents the average range residual over a one-day arc. Data affected by 
high plasma noise (intervals when SPE < 40°, indicated in gray) are not shown for clarity. The 
median and standard deviation of each time series are indicated in the legend. 
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Figure 10. Computed eccentricity for frozen orbits of varied inclination and semi-major 
axis (referenced to R = 2440 km to yield altitude h). Filled areas indicate that the frozen 
eccentricity is above the maximum acceptable eccentricity (dot-dashed line) and would lead to 
surface impact. Dashed lines indicate that the argument of pericenter cj is 90°, whereas solid lines 
indicate tu = 270° (pericenter over the southern hemisphere). Only the zonal gravity coefficients 
of HgM005 were considered. 
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Figure 11. Orbit evolution for the h = 500 km, i = 90° frozen orbit (found in Figure 10), 
as propagated over 8.5 years with the full HgM005 held and additional perturbations, (top) 
Periapsis altitude (in km), (bottom) Equinoctial element plot, showing the stability of the 
orbital elements. 
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Table 1. 

number ^ 2 . 


Model parameters used with ALMA [Spada, 2008], and the resulting tidal Love 


Elastic thickness 

Mantle rigidity 

FeS layer rigidity 

Computed k 2 

70 km 

75 GPa 

- 

0.46 

90 km 

75 GPa 

- 

0.46 

70 km 

65 GPa 

35 GPa 

0.55 

70 km 

95 GPa 

35 GPa 

0.52 

70 km 

10 GPa 

35 GPa 

0.62 
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Table 2. Kaula factor for Mercury inferred by scaling [Konopliv et ai, 2014] from other 
planetary bodies. 


Planetary 

body 

Kaula 

factor 

Mass®' 

(kg) 

Radius‘s 

(km) 

Scaled 

Kaula 

Reference 

Earth 

7.07 X 10“*^ 

5.97 X 10"^ 

6378.0 

4.96 X 10“3 

Lemoine et al. [1998] 

Moon 

3.6 X 10“4 

7.35 X 10^2 

1738.0 

6.94 X 10“5 

Lemoine et al. [2013] 

Venus 

1.2 X 10“5 

4.87 X 1024 

6051.0 

6.91 X 10“5 

Konopliv et al. [1999] 

Mars 

18.4 X 10“® 

6.42 X 1023 

3396.0 

1.86 X 10“4 

Lemoine et al. [2001] 

Vesta 

1.1 X 10“2 

2.59 X 1020 

265.0 

4.87 X 10“5 

Konopliv et al. [2014] 


^ From the JPL DE423 ephemeris [Folkner, 2010]. 

^ Reference radius for the gravity helds. Typically the same as the lAU convention, except 
for the Moon and Vesta. 
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Table 3 


Estimated values and formal uncertainties of selected low-degree coefficients 


Coefficient 

Value 

Formal uncertainty (cr) 

Ratio {C]^/a) 

C20 

-2.25045 X lO-*^ 

6.8 X 10-^*^ 

- 

C22 

1.24538 X 10“® 

4.4 X 10-10 

- 

C30 

-0.47659 X 10“5 

1.6 X 10“® 

- 

C21 

-1.61527 X 10“* 

3.8 X 10-10 

~ 42 

^21 

-1.36488 X 10“* 

3.9 X 10-10 

~ 35 

S22 

-2.09078 X 10“® 

2.3 X lO-o 

~ 9 
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Table 4. 


Brouwer-Clemence SET III parameter corrections, and associated ephemeris Keple 


rian element corrections. 


Parameter 

Unit 

Value 

Parameter 

Unit 

Value 

Aa/a 

- 

4.48 X 10“^^ 

Aa 

m 

0.26 

Ae 

- 

-8.83 X 10-1° 

Ae 

- 

-8.83 X 10-1° 

AMo + Auj 

O 

2.52 X 10-1° 

Ai 

O 

-1.40 X 10-° 

Ap 

o 

-5.27 X 10-1° 

An 

O 

2.30 X 10-° 

Aq 

o 

2.48 X 10-* 

Auj 

o 

-2.25 X 10-° 

eAuj 

o 

1.15 X 10-1° 

AM 

o 

- 
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